Aviram-Ratner rectifying mechanism for DNA base-pair sequencing through graphene 
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We demonstrate that biological molecules such as Watson-Crick DNA base pairs can behave as 
biological Aviram-Ratner electrical rectifiers because of the spatial separation and weak hydrogen 
bonding between the nucleobases. We have performed a parallel computational implementation of 
the ab-initio non-equilibrium Green's function (NEGF) theory to determine the electrical response 
of graphene — base-pair — graphene junctions. The results show an asymmetric (rectifying) current- 
voltage response for the Cytosine-Guanine base pair adsorbed on a graphene nanogap. In sharp 
contrast we find a symmetric response for the Thymine- Adenine case. We propose applying the 
asymmetry of the current-voltage response as a sensing criterion to the technological challenge of 
rapid DNA sequencing via graphene nanogaps. 
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I 

in 



c3 



i 

C 

o 
o 



(N 
> 

(N 
(N 
O 

o 

(N 

> 

X 

S-H 



I. INTRODUCTION 

Rapid and low-cost genome sequencing is one of the 
grand challenges of genome science today. The Sanger 
method [I] has served as the cornerstone for genome se- 
quence production since 1977, close to almost 30 years of 
tremendous utility. With the completion of the human 
genome sequence [H [5] , there is an imminent need in de- 
veloping new sequencing methodologies that will enable 
"personal genomics" or the routine study of our individ- 
ual genomes [H [5] . 

One potential candidate is nanopore sequencing [SHE], 
where a negatively-charged single-stranded (ss) DNA (in 
solution with counterions) is envisioned to translocate 
through a-hemolysin channels in lipid bilayers [ME] or 
through solid-state nanopores [9-14 , by a longitudinal 
electric field from one side of a membrane to the other. 
As the nucleotides of the DNA are migrated across the 
membrane, a significant fraction of ions are blocked from 
simultaneously entering the pore depending on their size. 
By continuously measuring the ionic blockade current, 
single molecules of DNA may be detected. Another pro- 
posed electronic approach is based on the intrinsic elec- 
tronic properties of the bases by, for instance, embedding 
electrodes within a nanopore to measure the transverse 
current through ssDNA as it translocates through the 
pore |14j . However, the construction of a nanopore with 
embedded electrodes remains a formidable challenge to 
the implementation and testing of this method. 

Recently, a systematic nanoelectrode-gated transverse 
electron-tunneling molecular detection concept with po- 
tential for rapid DNA sequencing has been proposed 
|15j . According to the nanoelectrode-gated molecular- 
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detection concept, it should be possible to obtain genetic 
sequence information by probing a DNA molecule base 
by base at a sub-nanometer scale. The nanoscale read- 
ing of DNA sequences is envisioned to take place at a 
nanogap defined by a pair of nanoelectrode tips as a DNA 
molecule moves through the gap base by base. The ra- 
tionale is that the four different nucleotide bases [ade- 
nine (A), thymine (T), guanine (G), and cytosine (C)] 
and their various sequences, each with a distinct chemi- 
cal composition and structure, should be associated with 
a specific signature of transverse tunneling current across 
the two tips. Theoretically, this new approach has the po- 
tential to sequence DNA at a maximal rate of 106 base 
pairs per second per detection gate. This method can 
be extended to parallel arrays of multiple nanoelectrode 
detection gates, thus magnifying the readout throughput 
by additional orders of magnitude, achieving estimated 
maximal rates of possibly hundred million (100 Mb) bases 
per second per device. Because the nanoelectrode-gate 
electron tunneling detection approach does not depend 
on the transient blockade of the ion current, one can use 
a wider range of detection gap sizes (1.5 - 5 nm) com- 
pared to the very small detection gap sizes (1.5 - 2.5 nm) 
required by the electrophysiology-like approaches. 

Despite the promise that transverse conduction mea- 
surements of DNA molecules holds for rapid sequencing, 
a comprehensive experimental study of transverse DNA 
conductance showing a variation of the conductance with 
base type has not been performed yet. Furthermore, the 
exact mechanism of electronic transport and its signature 
is debated [H HUffiS] . 

In this manuscript we study a novel approach for trans- 
verse charge transport through individual or a sequence 
of bases of double stranded (ds) DNA placed between 
the fringes of a graphene nanogap. The electronic sig- 
natures of the individual base or sequence of bases will 
display characteristic discrete peaks in the density of 
states, which will lead to a varying conductance finger- 
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print as the bias voltage is varied. The novel single-layer 
graphene-nanogap technique proposed here for transverse 
conductance measurements of molecules offers the fol- 
lowing desirable features compared to commonly-used 
techniques employing thin metallic (Au) wires: (1) ac- 
curate imaging of the location of the molecules inside the 
nanogap; (2) precise knowledge of the configuration of 
all atoms in contact with the molecule; (3) planar geom- 
etry which allows for the fabrication of parallel arrays of 
graphene nanogaps, thus magnifying the DNA sequenc- 
ing rate; and (4) tuning of the discrete molecular energy 
levels with respect to the leads' chemical potential. 

We demonstrate that the DNA base pair behaves as 
a biological Aviram-Ratner electrical rectifiers because 
of the spatial separation and weak bonding between the 
nucleobases. 



II. MATERIAL MODELS 
A. The two-level model 

In 1974 Aviram and Ratner [19] showed that a molecu- 
lar dimer composed of both an acceptor and a donor sub- 
units that are spatially separated by a sigma-bond bar- 
rier behave as a molecular rectifier. The two-level model 
in Figure la shows the relative position of the highest 
occupied HOMO (2, 4) and lowest unoccupied LUMO 
molecular orbitals (1, 3), with respect to the Fermi level 
Ep, for the acceptor and donor subunits of a generic 
Aviram-Ratner system. The application of a bias volt- 
age V across the junction causes a rigid shift of these 
energy levels by =F ¥ depending on whether they are spa- 
tially localized on the left or right side of the junction, 
as seen in Figure lb. The relative shift of the levels is 
responsible for the rectifying behavior, which is charac- 
terized by an asymmetric current-voltage (I — V) curve 
shown in Figure lc. 

For forward bias, the current rises around Vp = 
- min{Ai, A4} as either 1 or 4 enters to the Fermi energy 
window, highlighted in magenta. The schematic Figure 
lb shows a tunneling processes that is resonant through 
the donor level 4 but non-resonant across the acceptor 
subunit. For reverse bias, analogously, the current onset 
happens at Vr = imin{A 2 ,A 3 }. The donor-acceptor 
condition implies that the energies of the frontier levels 
of the subunits satisfy £2 < £4 < £1 < £3; consequently 
Vp < Vr, which causes the typical I — V asymmetry of 
Aviram-Ratner diodes. 



B. The single-base-pair model for transverse 
tunneling 

1. Longitudinal hopping vs. transverse tunneling 

The longitudinal electronic transport along DNA re- 
mains controversial and has not been conclusively deter- 



mined whether the DNA is metallic or insulating [2"Ull21| . 
Longitudinal transport directly involves the fluctuating 
chemical environment around the DNA backbones and 
thus can support multiple charge transfer mechanisms 
that arise from the small activation gaps induced by wa- 
ter and counterions. 

In contrast, the transverse electronic transport, per- 
pendicular to the dsDNA axis, is a simpler and less con- 
troversial process. It involves an insulating barrier (the 
hydrophobic core) and only a few discrete energy levels 
within the barrier, which belong to the base pair. 



2. The single-base-pair approximation 

The interaction between the stacked base pairs is neg- 
ligible (of the order of 0.01 eV for A-DNA and 0.1 eV for 
B-DNA [20 , 22 and therefore the total transverse current 
of DNA translocating through a nanopore can be well 
approximated as independent contributions from multi- 
ple channels. Each base pair temporarily located within 
the nanopore's electrodes (recognition region) constitutes 
an independent channel. In the case of zero-thickness 
graphene electrodes, our model approximates transverse 
tunneling as through a single base pair that is decou- 
pled from its neighbors, as shown in Figure 2. The ionic 
environment around the backbone (ions, counterions, sol- 
vent) and the dynamics of the translocating process are 
necessary for a complete description of in vivo DNA; 
nonetheless, they involve computationally intense calcu- 
lations which are out of reach currently employing solely 
ab initio calculations. In addition, as discussed later, 
both experiment and theoretical calculations have shown 
that the transverse transport primarily depends on the 
nature of the nucleobases rather than on the in vivo en- 
vironment. Our proof-of-concept model attempts to ad- 
dress the underlying physics of transverse transport and 
focuses on single base pairs. 

3. Effect of the backbone, solvent, and counterions 

Stability of the energy levels of the base pairs: 
Regarding the solvent, fluctuations of the surrounding 
water molecules are known to have little effect on the 
transverse current for the case of ssDNA, amounting to 
small modulation of its magnitude [18| . More impor- 
tantly, the energy levels for dsDNA are expected to be 
more stable against external perturbations compared to 
those of the ssDNA each base pair is protected in- 

side the hydrophobic core that is further stabilized by the 
interaction between backbones and counterions. In cases 
when water can enter into the DNA structure through 
damaged sites, it can induce small activation gaps around 
the Fermi level [23] . 

In regard to the backbone, the discrete energy levels 
of the base-pairs, relevant to transport, are electronically 
isolated from the continuous spectrum of the backbone's 
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phosphate and sugar (i.e. off-resonance 1.24]), which in 
turn have a much larger energy gap of ~ 6.5 eV |25) and 
hence do not hybridize significantly with the base pair 
within the relevant energy window |26l 12?j . The presence 
of countcrions (e.g. Na+, Mg+) is more critical, how- 
ever, since they introduce energy levels in the base-pair's 
HOMO-LUMO gap by hole doping [201 [2SJ. Nonetheless, 
these states are localized outside the backbones and could 
induce electron hopping (in conjunction with nearby wa- 
ter) along the longitudinal direction but are unlikely to 
contribute to transverse tunneling, which proceeds pri- 
marily through the hydrogen-bonded base pairs. Recent 
experiments find that the leakage current due to these lo- 
calized states amounts to a small offset and noise in the 
measured tunneling current and is readily separable from 
the main contribution coming from tunneling through the 
nucleobases [29] . 

Furthermore, shift of the discrete levels due to geomet- 
rical deviations from the optimal adsorption (7r-stacking) 
of the base pair on the graphene electrodes caused by 
thermal fluctuations or backbone-strain relaxation (not 
considering translocation) is also negligible. Combined 
molecular dynamics (MD) and density functional theory 
(DFT) statistical analysis of nucleobase adsorption on 
graphene [26 , 30J and CNT [3T] confirm that this type of 
perturbation has small effect and that the relative posi- 
tions of the nucleobase's levels are locked, almost inde- 
pendently from the orientation between the nucleobases 
and the electrode. 

Alignment between the energy levels and the electrode's 
Fermi level: 

Another relevant factor in the Aviram-Ratner mecha- 
nism, besides the stability of the discrete energy levels 
discussed above, is the alignment between these levels 
and the Fermi level of the electrodes Ep. Such align- 
ment determines the parameters Aj (-1=1-4) and is ex- 
pected to be susceptible to the relative geometrical orien- 
tation between the hydrophobic core and the electrodes. 
At the equilibrium optimal adsorption, the alignment 
is mainly determined by a van der Waals interaction. 
However, far from equilibrium the alignment is gener- 
ally determined by the long-range Coulomb interaction 
that couples the nucleobase to the electrodes. The dy- 
namical displacements during DNA translocation cre- 
ates a random electrostatic environment that overlaps 
with the weaker nucleobase-electrode interaction which 
in turn can cause dynamical fluctuations in the HOMO- 
EF alignment and thus in Aj. Since the nucleobase- 
electrode interaction is ultimately dependent on the elec- 
trical polarizability of each individual nucleobases (G > 
A > T > C) (32H33], the values of A t should also fol- 
low a consistent trend as long as the polarizability of 
the nucleobases can be consistently detected in the midst 
of the environmental electronic noise. Although appar- 
ently counterintuitive, MD simulations have proven that 
such weak nucleobase-electrode interaction can be distin- 
guished from the chaotic environmental electronic noise 
[35j . Statistical averaging or increase of residence time 



of the nucleotide inside the recognition region is needed 
to manage the electronic noise |36j . Distinguishing such 
a weak interaction is equivalent to measuring the capac- 
itance of each nucleobase, which is another promising 
technique for fast DNA sequencing with the single-base 
resolution [37] . 

Experimental reproducible determination of the energy 
levels: 

The HOMO-Ep alignment is primarily pinned by the 
intrinsic electronic signature of each nucleobase. De- 
spite the expected random electrostatic environment of in 
vivo DNA, electrostatic cancelation effects arc expected 
[26], resulting in a more electrostatically homogeneous 
medium, and consequently a less fluctuating HOMO-£i? 
alignment. In that regard, by measuring the transverse- 
tunneling current (with scanning tunneling spectroscopy 
STS), experiments have successfully pinpointed the value 
of the HOMO level with respect to Ep within a 10% and 
4% standard deviation for poly(CG) and poly(TA) ds- 
DNA, respectively [551150], 

Moreover, the same STS technique has also repro- 
ducibly measured the nucleotide's HOMO levels in ss- 
DNA [3H HO] , despite the fact that ssDNA is conforma- 
tional less stable than dsDNA. Furthermore, a more re- 
cent experiment using a simple break junction instead 
of a STM has further proved that such identification is 
possible, a step closer to our proposed junction model 

In short, recent STS measurements have shown that it 
is possible to pin the energy levels of nucleobases with 
encouraging degree of reproducibility, despite the pres- 
ence of backbone, solvent, and counterions. Moreover, 
an Aviram-Ratner rectifying shape of the transverse cur- 
rent, similar to that predicted in this manuscript, has al- 
ready been experimentally observed in poly(CG) dsDNA 
[25]. 



III. COMPUTATIONAL METHODS 

All the electronic structure calculations are performed 
at the ab initio level using the SIESTA method [42] 
with the generalized-gradient Perdew-Burke-Ernzerhof 
exchange-correlation density functional [33] . The valence 
electrons are expanded using localized basis of double-£ 
with polarization orbitals (DZP). We found that polar- 
ization orbitals are especially needed to screen the ap- 
plied bias voltage along the electrodes. Troullier-Martin 
norm-conserving pseudopotentials [44j in the Kleinman- 
Bylander form |45| are used for the core electrons. 

The electronic density is considered converged after the 
maximum difference between the density matrices of two 
consecutive cycles is smaller than 10~ 5 . 

The positions of the atoms in the unit cell of the 
nanoribbon (a = 4.92A), the reconstruction of the 
nanogap, and the adsorption geometry of the base-pairs 
on the nanogap were self-consistently optimized by con- 
jugate gradient until the forces are less than 0.01 eV/A. 
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For the isolated base pairs, ample vacuum separation 
was included inside the unit cells to avoid interaction 
with mirror images. In the case of the infinite nanorib- 
bon, the calculations used 80 k-points along the periodic 
direction. For the case of the finite C-G and T-A junc- 
tions, in order to speed up the convergence of the wave- 
functions, the boundary atoms in the unit cell are made 
to match periodically along the longitudinal direction; 
ample vacuum is included in the other two directions. 



For the case of an "infinite" junction, the retarded 
Green's function is obtained by including the effect of 
the left (L) and right (R) semi-infinite electrodes through 
self-energy corrections Er^m- to the the finite system H 
|46| . resulting in the following expression 

G+ = G+S = [e+S -H- S L (e+) - E H (e+)]- 1 5 (7) 
This leads to the following expression 



A. Non-equilibrium Green's function (NEGF) 

The central part of the NEGF method involves calcu- 
lating the electron density of an infinite junction under 
bias voltage. This subsection presents a self-contained 
derivation leading into a matrix representation of the 
non-equilibrium electron density, given in Eq. |12| A 
matrix representation is important since all necessary 
quantum-mechanical quantities such as Hamiltonian and 
overlap operators can be readily extracted from any 
localized-orbital DFT package in matrix form. 

The spectral-function operator A of a system with 
Hamiltonian operator H is defined as 



A{e) = 2irS(e - H) 



(1) 



From the Sokhotsky-Weierstrass theorem, the retarded 
(+) and advanced (-) Green's function operators G ± (e) 
can be decomposed as: 



G ± (e) = lim 



i-yo e — H ±rji 



V- 



1 



s T«4-5), (2) 



H 



where V represents the Cauchy principal- value of an asso- 
ciated integral. Using this relation, the spectral-function 
operator reduces to 



A(e) ee 2tt<S(£ -H) = i{G + - G' 



(3) 



(G+) 1 -{G~ 



r 1 =e + S-H-Z L (e + )-i; R (e + ) 
-(£-£-#-££(£-)-£«(£-)) 
= (e+ -e-)S + i{T L +T R ) 



where the definition of the broadening matrix r^ m = 
i(T,r L R }(e + ) — ^{l,r}( £ ~)) was used in the last line . 
Multiplying both sides of the equation through the left 
and right by —iG + and G _ , respectively, one finds 



i(G+ - G~) = 2r}G + SG~ + G+(T L + r R )G~ 



(8) 



The Green's function matrix in an orthogonal basis is 
obtained through G = GS [52 . Combining Eq. [3] and 
Eq. [8j one finds 



A = i(G + -G~) (9a) 
A = 2 V G + SG-S + G + {T L +r R )G-S (9b) 

By definition, the density of states of the junction per 
spin channel is given by 



DOS = Tr( 



Tr( 



AS~ 
2tt 



(10a) 



In the absence of coupling of the central region to the 
electrodes (r^ = F R = 0), it follows from Eq. 9b that 
the density of states of the junction reduces to |47j 



The density-matrix operator D is defined as the Fermi 
function f(H-f-i) = [exp ((H - fi)/(k B T))} + 1]" 1 of the 
Hamiltonian operator H at the "equilibrium" chemical 
potential /i 



D{n) = i[H-n) (4) 
f(e - (i)S(e - H)de (5) 



1 

27 



+ 00 



def(E-n)A(s) 



(6) 



The bars in the last expression denote matrices ex- 
panded in an orthogonal basis set. Otherwise, all matri- 
ces are given non-orthogonal basis by default. 



DOS 



- rf 



-An^ n ) (11) 



where c„ = tp^S/fin- £n and <j> n are the generalized eigen- 
values and eigenvectors of the finite system defined by H 
and S. From Eq. [Tl] it can be verified that the density 
of states of a uncoupled junction reduces to a series of 
discrete energy levels represented by Lorentzian distri- 
butions with 77 — ¥ 0. 

In the presence of coupling, the energy levels of the 
central region are no longer discrete but continuous in 
energy, as schematically depicted by the broad horizon- 
tal lines in Figure From Eq. |9b[ the "scattered" 
electronic states of the infinite junction fall into one of 
the following cases: (i) states that couple the central 
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region and left lead (the term G + TlG~ S); (ii) the cen- 
tral region and the right lead (the term G + TrG~ S)\ the 
states that couple the central region and both leads are 
projected into the previous two cases; (Hi) states local- 
ized in the central region that do not couple to any of the 
leads (the term 2r\G + SG~ S); and (iv) states that couple 
both leads but not the central region, which are notice- 
ably absent in Eq. |9b| as a consequence of the underlying 
assumption of non-iteracting electrodes. 

The density matrix gives the electronic charge of the 
junction and is found by filling (integrating) the scattered 
states over the energy range where reservoir electrons 
are available in the electrodes, as expressed by Eq. [6] 
Without applied bias voltage, the states corresponding 
to cases i and ii are evenly filled up to the chemical 
potential of the junction (/i = Ep). 

Under applied V, however, states i and ii are unevenly 
filled since they are required to reach local equilibrium 
with the chemical potential of the left (/j,l = Ep — ^) or 
right (fi R = E F + ^) electrode, respectively. As schemat- 
ically shown in Figure [3^,, the total electronic charge of 
the central region derives from the states in the green 
and gray areas . Therefore, the charge is said to be 
decomposed into an "equilibrium" (green) and a "non- 
equilibrium" (gray) contribution. 

The equilibrium contribution D eq is analogous to a fic- 
titious zero-bias problem (given by Eq. [6]) where both 
electrodes have the same chemical potential Then, 
D eq — D(/il). In this fictitious zero-bias problem some 
information of the applied bias voltage is implicitly in- 
cluded in the shifted electronic structure of the right 
electrode. The non-equilibrium contribution D non - eq re- 
sults from all the states inside the gray area, which do 
not couple to the left electrode and therefore can be ob- 
tained from the spectral function in Eq. |9b|when Tl = 0. 
Following Eq. [6j we compute this contribution by inte- 
grating A| ri _ over the f(e - fi R ) - f(e - fi L ) energy 
window. Combining Eq. [6] and Eq. |9b[ the density ma- 
trix for the "infinite" junction under bias voltage is: 



D(h<l,Hb,) — D eq + D non - eq 

-t-oo 



1 

2^ 



dei{G+ -G-)f{e- ^l) 



+ 00 

^ J de[2r]G + SG-S + G + r R G-S} 



(12) 
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FIG. 1: Two-level model of an Aviram-Ratner molecular rec- 
tifier, (a) 2 and 1 (4 and 3) are the HOMO and LUMO levels 
of the acceptor (donor) subunit. (b) Rigid shift of the energy 
levels under an external bias voltage V = \&-4, applied across 
the electrode's reservoir states (shown in blue). Both, the in- 
terfacial and the internal energy barriers are shown in gray. 
The Fermi energy window due to V is shown in magenta, 
(c) Asymmetric current-voltage curve of an Aviram-Ratner 
diode; the forward (Vf) and reverse (Vr) threshold voltages 
for the onset of the current depend directly on the relative po- 
sition of the frontier energy levels, (d) Ab-initio energy levels 
for the isolated C-G and T-A base pairs, obtained from the 
projection of the density of states (PDOS) on the pyrimidine 
(blue) and purine (red) subunits. 



IV. RESULTS AND DISCUSSION 

A. Isolated Cytosine-Guanine (C-G) and 
Thymine- Adenine (T-A) base pairs 

The original Aviram-Ratner rectifier employs a sigma- 
bond internal barrier between the 7r-systems of the sub- 
units. For the case of DNA base pairs, both subunits 
are weakly coupled by hydrogen bonds that act as in- 
ternal insulating barriers. Figure Id shows the density 
of states of isolated C-G and T-A molecules projected 
on their subunits. The isolated molecules arc calculated 
without electrodes but using the geometry they adopt in 
the junction environment. 

First, the presence of sharp peaks indicates minimal 
hybridization between the subunits, which is a direct con- 
sequence of the weak hydrogen-bond coupling. There- 
fore, the electronic structure of the base pair can be ap- 
proximately viewed as the superposition of the two sub- 
units, which fulfills the barrier condition. Second, the or- 
dering of the computed energy levels (£2 < £4 < £i < £3) 
corresponds to that of an acceptor-donor pair. There- 
fore, both C-G and T-A base pairs are biological Aviram- 
Ratner rectifiers. Our ab-initio values of Aj (z=l-4) for 
the isolated C-G and A-T are 1.61 eV, 2.11 eV, 3.27 eV, 
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FIG. 2: DNA base-pair junctions. Top and side view of the 
C-G (a) and T-A (c) junctions. The density of states are pro- 
jected on the purine (red) and pyrimidine (blue) components 
of the C-G (b) and T-A (d) junctions. The projected eigen- 
states are shown for reference on the middle panels of (b) and 
(d); the colored levels correspond to the four relevant ones in 
the two-level model. The relevant projected eigenstate 4* is 
located at -0.38(-0.92) eV for the C-G (T-A) junction. The 
unit cell and bandstructure for the reconstructed ribbon are 
shown in (e) and (f); the ribbon exhibits a forbidden-energy 
region, below the Fermi level Ef = eV, highlighted in green. 



0.41 eV, and 1.79 eV, 1.76 eV, 2.16 eV, 1.6 eV, respec- 
tively. For both base pairs A4 < Ai and A2 < A3. 
Consequently, the forward (Vp = j^i) and reverse 
(Vr = -A 2 ) threshold biases are determined primarily 
by the HOMOs of the purine (level 4) and pyrimidine 
(level 2) subunits, respectively. Moreover, because of the 
large difference between the C and G HOMO energies 
(A 2 A 4 ), C-G is expected to exhibit a noticeably 
more asymmetric I — V curve (pP = 0.19) than T-A 

(v = 0-91)- We propose that this difference in I — V 
asymmetries can be used as an electrical probe for dis- 
cerning between C-G and T-A pair. 



B. GNR — base-pair — GNR junctions 

The presence of metallic electrodes in the molecular 
junction plays an important role in modifying the en- 
ergy levels of the base pair through hybridization; it 
may result in degradation or even loss of the rectifying 
mechanism. Because of the low density of states at Ep, 
carbon-based electrodes have been proposed as highly 
sensitive detectors with molecular resolution suitable for 
rapid DNA sequencing [26J and as non-intrusive intercon- 



necting wires for electronics 48 . However, low density of 
states at Ep also indicates electrodes with rapidly decay- 
ing wavefunctions and therefore small variations to the 
adsorption geometry of the base pair, while translocat- 
ing through the nanogap, can induce large fluctuations of 
the tunneling current; nonetheless, nonlinear techniques 
have been proposed to overcome that difficulty |4"9] . 

We use reconstructed graphene nanoribbons (GNRs) 
as electrodes. The reconstructed edges are composed of 
heptagons and pentagons and increase the dispersion of 
the otherwise flat low-energy bands; therefore, it elimi- 
nates the well-known spin-induced bandgap of unrecon- 
structed zigzag ribbons, thus rendering the reconstructed 
GNR metallic. Reconstructed ribbons have been pre- 
dicted to be energetically favorable [SO]- We construct 
a nanogap of ~ 5A in the reconstructed GNR, whose C 
edge atoms are passivated with H. The geometry of the 
nanogap is further relaxed and kept fixed in subsequent 
DNA adsorption, as seen in Figure 2a. Our transport 
calculations confirm that the tunneling current through 
the bare nanogap is negligible. 

The C-G and T-A base pairs were allowed to relax on 
top of the graphene nanogap. In both cases, the pyrimi- 
dine subunit (G or A) adsorbs almost parallel above the 
plane of the ribbon at ~ 3A following the same Bernal's 
AB stacking of graphite that minimizes 7r — ir repulsion 
and indicates physisorption via electrostatic coupling, as 
previously reported for nucleobases [3SJ [5T] . In the ab- 
sence of bias, the electronic structure of the "infinite" 
GNR — base-pair — GNR junction is obtained following 
the Green's function formalism, where the left and right 
semi-infinite electrodes are "mathematically attached" to 
a finite "central region" \$<o\ . Figure [2] shows the central 
regions for the resulting C-G and T-A junctions. Because 
of the short screening length of the carbon electrodes, the 
central region includes multiple screening layers and has 
a length equivalent to 15 GNR unit-cells, as shown in 
Figure [4£i. 

In Figure 2b we show the projected electronic density 
of states (PDOS) of the infinite junction on G (red) and 
on C (blue) [52]. The eigenenergies of the junction "pro- 
jected" on the C-G base pair are shown with vertical 
lines. First, it is observed that the projected eigenener- 
gies keep the same ordering than the eigenenergies of the 
isolated C-G (Figure Id). Second, the PDOS consists of 
a sequence of narrow peaks that closely follow the pro- 
jected eigenstates, which confirms that the electrodes do 
not significantly alter the molecular properties of the iso- 
lated base pair because of weak electrostatic coupling. 
Both observations further indicate that the junction may 
retain the Aviram-Ratner diode behavior predicted for 
an isolated base pair. 

The Lorentzian-likc G HOMO peak, around Ep = 
eV, clearly derives from 4*, which corresponds to 4 of 
the isolated base pair (asterisks denote a junction's eigen- 
state "projected" on the base pair). Therefore, according 
to the two-level model, the HOMO peak is expected to 
dominate the rectifying properties of the junction. The 
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PDOS on C is negligible, indicating that the HOMO peak 
is spatially confined to G, analogous to the isolated-C-G 
case. 



is closed with a counterclockwise arc C in the upper- 
half complex plane extending from +00 to —00. From 
Cauchy's residue theorem applied to the integral of I 
along the described closed contour R + C, we have 



C. Confirmation of the Aviram-Ratner mechanism 

For most cases, it suffices to calculate the electri- 
cal response of a molecular junction only in the linear- 
response regime (Landauer-Buttiker) , where the trans- 
mission channels are considered independent of the ap- 
plied bias voltage. However, the frontier transmission 
channels (eigenstates) in Aviram-Ratner rectifiers inher- 
ently undergo electronic changes with the applied V and, 
therefore, it requires an explicit recalculation of the elec- 
tronic charge inside the junction at every bias-voltage 
point (non-linear or Keldysh regime). Following our pre- 
vious developments |46j[52], we have implemented an in- 
house parallel software for the computation of the elec- 
tronic transport properties under non-equilibrium condi- 
tions. 

The electronic charge of the junction is represented in 
matrix form by the density matrix D. The density ma- 
trix for the central region connected to two semi-infinite 
electrodes and under applied bias voltage V was derived 
in Eq. [12] In a non-orthogonal basis-set it is expressed 
as 



D{hl, Mil) — D eq + D non _ eq 
+00 

de(G + -G-)f(e-fi L ) 



+00 



1 

2^T 



4-oo unbound 



bound 



(13) 



-!- / de[2r)G + SG- + G + T R G-] 
2n J 



x [/(£- m)~ /(e- Mi)], 

The Hamiltonian matrix H of the central region un- 
der applied V is computed using a standard DFT solver 
(SIESTA 2.0.2 42J), modified to introduce V as an exter- 
nal Hartree potential. Starting from an initial guess, the 
density matrix is calculated according to Eq. [13] and fed 
back to the DFT solver, which yields an updated Hamil- 
tonian matrix. The process is repeated iteratively until 
self-consistent convergence is achieved. 



D(H) i 



DFTsolv 



NEGF 



H{D) 



(14) 



The equilibrium contribution D eq in Eq. 13 is deter- 
mined by integrating the states inside the green area 
schematically shown in Figure [3^,. Because of the an- 
alytic properties of the Green's function, the integral of 
1(e) = (G + (e) — G~(e))f(e — hl) can be analytically 
continued to the complex plane, away from the real-axis 
singularities. The initial real-axis path R — [—00; +00] 



dzX(z) 



R+c 



del(s)+ / dzL(z 



2iri Res(I, z v ), 



(15) 

where Res(I,z v ) is the residue of I(z) evaluated at the 
singularities z v inside the closed contour. Those singular- 
ities are assumed to be first order and to come only from 
the poles of the Fermi function, z v = [iL+i(2is + l)Trk B T. 
The residues of the Fermi function are Res(f, z v ) = 
-k B T. 

Res(T, z u ) = lim (z — z l/ )X(z) 

= lim (z-z v )f(z-fi L )(G+(z)-G-(z)) 
= [ lim (z - z„)f{z - fi L )](G+(z u ) - G~(z v )) 

Z—¥Z U 

= Res(f, z v )(G + (z u )-G-(z u )) 
= -k B T(G+{z y )-G-(z u )), 

(16) 

In practice the counterclockwise path C is approxi- 
mated to C = C\ in Figure [3b. The equilibrium contri- 



bution D eq in Eq. 13 reduces to 



D P 



dzl(z) 



k B Tj2lG + (z»)-G-(z„)}. 



27r JC, 

(17) 

The path d is defined by 7 = 20k B T, S = 10fc B T, 
A = 2irn p k B T, and E min — /i^ — 30 eV, where n p — 5 
is the number of Fermi-function poles considered, k B 
is the Boltzmann constant, and T = 300A". The in- 
tegrals along the line and arc of C\ are performed us- 
ing a 10-point trapezoidal integration and a 85-point 
Gaussian-Legendre quadrature, respectively. The in- 
finitesimal i] = 10~ 8 eV. The non-equilibrium contribu- 
tion D non ^ eq in Eq. 13 involves an integration along the 



path C2, infinitcsimaliy away from the real-energy axis, 
from {v> R + 7k B T)+i£ to (fi L -7k B T)+i£, with£ = 10~ 8 
eV. 

Having determined a converged electron density, 
jjconvg ^ £ Qr foe biased junction, its electronic proper- 
ties such as PDOS, transmission-probability function 
(TF(e, V)), and current-voltage characteristics (I — V) 
can be determined from H(D conv9 ) and the well-known 
linear-response transport equations. 



TF(e,V)=Tr(T L G+T R G-), 



(18) 



and 



I(V) 



J deTF(e,V)\f(e-ii R )-f{e-fi L )] (19) 



6eui| 

FIG. 3: (a) Schematic representation of the electronic states 
of an infinite junction under an applied external bias voltage 
V. The finite broadening of the discrete energy levels of the 
central region (horizontal lines) schematically accounts for the 
coupling to both semi-infinite electrodes. Under bias, the 
reservoirs of electronic states of the left and right electrodes 
(blue regions) are assumed to shift rigidly around the Fermi 
level (Ef). (b) The complex-energy-plane paths Ci and C2 
used for the integration of the states in the green and gray 
area, respectively. The arc in Ci is centered on the real-energy 
axis. The blue dots represent the poles of the Fermi function 
of the left electrode that are enclosed by C\ and the real axis. 



where n a = 2 accounts for the spin degeneracy while all 
remaining terms are defined in the Computational Meth- 
ods section. 

Figure [4^, shows the drop of the electrostatic poten- 
tial for V = —0.5 V along the C-G junction. Note that 
the voltage drop occurs mainly inside the nanogap, con- 
firming that the external potential is well screened at 
the boundaries of the central region. Moreover, as seen 
in Figure |4j:, the electronic states of the right electrode 
shift in energy with applied voltage at a rate of ~ 0.44 
eV/V, which is only slightly below the rigid-shift rate of 
0.5 eV/V and thus corroborates negligible voltage drop 
inside the right lead. Because of electrostatic coupling to 
the electrodes, both subunits of the base pair follow the 
local electrostatic potential of the GNR on which they 
are adsorbed. Figure |4Jd shows the G HOMO peak shift- 
ing in concert with applied bias, although at a slightly 
lower rate (~ 0.36 eV/V) than the right electrode. The 
rate mismatch indicates a small voltage drop building up 
at the interface between G and the tip of the electrode. 

According to the two-level model, the onset of the 
forward current for the C-G junction is expected at 
V F = ~ 0.38 V as the HOMO peak enters the 

tunneling energy window. Accordingly, the transport cal- 
culations find the onset at ~ 0.35 V, as seen in the blue 
curve in Figure [5]d. Under reverse bias, conversely, the 
HOMO peak shifts away from the tunneling window and 
consequently, since level 1* does not participate, the re- 
sulting current is low and due only to nonresonant "leak- 
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FIG. 4: (a) Bias-voltage drop along the central region of 
the C-G junction. The left and right leads in our model con- 
tains 6 and 7 GNR unit-cells, respectively. The applied bias 
V = —0.5 V drops entirely around the nanogap and not along 
the electrodes, (b) Bias evolution of the density of states 
projected on Guanine (red). The broad peak, which shifts 
with bias, corresponds to Guanine's HOMO. Fermi windows 
(magenta) at every bias voltage are shown as reference, (c) 
Bias-voltage evolution of the density of states projected on 
the right electrode. 

age" tunneling (see later discussion on reverse-bias leak- 
age current). As a result, the I — V curve of the C-G 
junction adopts the rectifying shape for a diode. It is im- 
portant to emphasize that this I — V feature can be used 
to discern whether a C-G or a G-C base pair is translo- 
cating in the GNR nanogap. Similar differentiation has 
been reported in base pairs chemically adsorbed on gold 
electrodes [2"5] . 

D. Behavior of the T-A junction 

Analogous to C-G and according to the two-level 
model, the T-A junction should also exhibit a rectify- 
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ing I—V behavior dominated by the HOMO level 4* of 
A. However, for the T-A case, 4* is well below the Fermi 
energy in Figure 2d and hence could only be reached at 
rather high biases (Vf — -A4* ~ 0.92 V). Accordingly, 
the electrical conductance of the T-A junction is found 
to be much lower than C-G for the forward-bias points 
shown in Figure [5j which confirms that the Adenine's 
resonant level 4* is not contributing to conduction. 



1. Leakage current for reverse and low biases 

The reverse-bias current for the T-A is small and com- 
parable to C-G. Such "leakage" current arises from dou- 
ble non-resonant tunneling through both subunits and is 
analogous to Bardeen's tunneling theory for a scanning 
tunneling microscope tip (Adenine) on top of a surface 
(Thymine) 53J. The leakage tunneling in the T-A junc- 
tion is determined by the number of electronic states 
available on both T and A. Moreover, the number of 
states on T is ^ 2 two orders of magnitude lower than 
those on A; thus, T constitutes the limiting factor in the 
leakage-tunneling process. Bardeen's tunneling mecha- 
nism offers a simple explanation for the sigmoidal shape 
of the T-A current based only on the mesa-like shape of 
the zero-bias T PDOS . At low energies (\E - E F \ < 0.1 
eV), the mesa- like T PDOS is roughly constant (Figure 
[5^,) which , upon integration, yields the linear (ohmic) re- 
gion seen in Figure [5]d- red for |V| < 0.2 V. Moreover, in 
contrast to A, the T PDOS are spatially localized on the 
left side of the nanogap and hence shifts towards lower 
energies under applied bias V, as seen in Figure 5a. The 
number of states (area within the Fermi energy window) 
progressively increases with bias (linear region), until it 
saturates at V ~ +0.3 V and —0.7 V, as seen in Fig- 
ure [5^,, marking the forward- and reverse-bias saturation 
points of the I — V red curve. 

As mentioned before, for V > Vp ~ 0.92 V, one may 
expect larger resonant-tunneling currents mediated by 
level 4* of A. However, there are two mechanisms that 
may prevent such trend: 



2. Inhibited hybridization of the HOMO level 

In contrast to G, the level 4* of A does not fully hy- 
bridize with the states of the right electrode, since 4* falls 
inside the forbidden energy region of the GNR (green 
highlighted area in Figure 2f). The lack of hybridiza- 
tion between 4* and the right electrode hinders resonant 
tunneling and dramatically reduces the electrical current 
through the junction. Nonetheless, because of the close 
proximity of 4* to the GNR states above the forbidden 
region (i.e. distance to the upper boundary of the for- 
bidden region), partial hybridization is still possible at 
those energies. The partial or inhibited hybridization 
yields a truncated Lorentian-like HOMO peak in Fig- 
ure [2k, which is centered at 4*. Accounting for a much 



weaker hybridization, the truncated peak in Figure 2d is 
~ 8 times smaller than the HOMO peak in Figure 2b. 
Finally, a small surge in current may still be expected as 
the truncated peak enters the Fermi tunneling window, 
at large biases, analogous to the C-G HOMO peak in Fig- 
ure]^; however, this trend is further prevented because 
the truncated peak progressively decays with increasing 
bias voltage, as seen in Figure [6] 

3. Voltage drop at the Adenine/ electrode interface 

In order to elucidate the origin of the anomalous be- 
havior of the A HOMO peak, we followed the evolution 
of the peak and the level 4* with bias voltage. As seen 
in Figure 6, both 4* and the HOMO peak shift toward 
higher energies, but at different rates. The sharp peak 
corresponding to 4* progressively lags in relation to the 
green forbidden region. This indicates a small voltage 
drop building up at the interface between A and the tip 
of the right electrode. With increasing bias, 4* lags to- 
wards the lower boundary of the forbidden region. The 
increasing mismatch between 4* and the allowed states 
(denoted by the growing arrows in the figure) further 
prevents hybridization resulting in the observed decay of 
the A HOMO peak. Consequently, the expected surge of 
current at ~ Vf is considerably averted. 

It is important to emphasize that the unexpected de- 
parture of the T-A junction from the Aviram-Ratner be- 
havior is solely due to the fortuitous location of Ade- 
nine's level 4* within the forbidden energy bandgap. 
The bandgap itself is an artifact of the strong quan- 
tum confinement due to the ultra narrow graphene elec- 
trodes (width of ~ 1.2 nm) employed for computational 
tractability. The quantum confinement in wider ribbons 
is weaker and the dispersion of the low-energy bands in- 
creases, effectively eliminating such bandgap 50 . 

V. CONCLUSIONS 

We have demonstrated that in the presence of GNR 
electrodes, Watson-Crick DNA base pairs still act as bi- 
ological Aviram-Ratner rectifiers because of the spatial 
separation and weak bonding between the two nucle- 
obases. 

The C-G pair junction shows a noticeable rectifying be- 
havior with highly asymmetric I—V curve, in close agree- 
ment with the Aviram-Ratner mechanism for molecular 
rectifiers. We propose that the I — V asymmetry can be 
used as an electronic handle to discern C-G from G-C in 
the context of DNA base-pair sequencing through trans- 
verse electrical measurements using graphene nanopores. 

In sharp contrast, we find that the T-A junction ex- 
hibits a rather symmetric I—V with sigmoidal shape, 
at least for moderate voltages |V| < 0.9 V. The electri- 
cal characteristics of both C-G and T-A junctions are 
dominated by the HOMO level of the purine nucleobase 
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(either Guanine or Thymine). Because of the markedly 
different relative positions of their HOMO peaks with 
respect to E F , the C-G junction exhibits larger current 
than the T-A junction for positive bias. For reverse bias, 
however, both junctions show "leakage" currents of com- 
parable magnitude. 
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FIG. 5: (a) Evolution of the density of states projected on 
Thymine for various bias, (b) Current- Voltage (7 — V) curves 
for the C-G (blue) and T-A (red) junctions. 




FIG. 6: Evolution of the truncated Adenine HOMO peak and 
the peak corresponding to 4*, both shown in Figure|2ji, under 
an applied bias voltage V . The arrows show the progressive 
energy mismatch (0.09 eV, 0.13 eV, 0.16 eV, and 0.18 eV) 
between both peaks. For visual aid, the forbidden energy 
region of the right electrode (green) and the zero-bias PDOS 
(dotted lines) are shown in each panel under a rigid energy 
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